Spatial dynamics of malaria transmission

The Ross-Macdonald model has exerted enormous influence over the study of malaria transmission dynamics and control, but it lacked features to describe parasite dispersal, travel, and other important aspects of heterogeneous transmission. Here, we present a patch-based differential equation modeling framework that extends the Ross-Macdonald model with sufficient skill and complexity to support planning, monitoring and evaluation for Plasmodium falciparum malaria control. We designed a generic interface for building structured, spatial models of malaria transmission based on a new algorithm for mosquito blood feeding. We developed new algorithms to simulate adult mosquito demography, dispersal, and egg laying in response to resource availability. The core dynamical components describing mosquito ecology and malaria transmission were decomposed, redesigned and reassembled into a modular framework. Structural elements in the framework—human population strata, patches, and aquatic habitats—interact through a flexible design that facilitates construction of ensembles of models with scalable complexity to support robust analytics for malaria policy and adaptive malaria control. We propose updated definitions for the human biting rate and entomological inoculation rates. We present new formulas to describe parasite dispersal and spatial dynamics under steady state conditions, including the human biting rates, parasite dispersal, the “vectorial capacity matrix,” a human transmitting capacity distribution matrix, and threshold conditions. An R package that implements the framework, solves the differential equations, and computes spatial metrics for models developed in this framework has been developed. Development of the model and metrics have focused on malaria, but since the framework is modular, the same ideas and software can be applied to other mosquito-borne pathogen systems.


Reproductive Numbers
Here we summarize the notation used in the text to define reproductive numbers. For a longer discussion, see Supplement 8.
R 0 -The Basic Reproductive Number is a scalar. Here, we have defined R 0 as the spectral radius of one of the parasite's next generation matrices, either R or Z. It describes transmission in a population that has not been modified by immunity or control, and it counts every infectious bite that would cause an infection, even if the host is already infected.
R C -The Adjusted Reproductive Number is a scalar. In fact, R C is like R 0 in every way, except one -R C has been modified by malaria control. We think of R C as describing a family of numbers, of which R 0 is a special case where C = 0. The difference between R 0 and R C are related to the interpretation of D and V.
To compute R C , D and D should be computed for a population that is lacking immunity, which may require changing some of the parameters describing transmission from strata that are immune.
R E -The Endemic Reproductive Number is a scalar. It describes transmission in a population that has been modified by immunity and control, but like R C , it counts every infectious bite that would cause an infection, even if the host is already infected. The endemic number thus includes effect sizes of immunity on transmission, and R C /R E .
R e -The Effective Reproductive Number is a scalar. It describes transmission in a population that has been modified by immunity and control, but it does not count redundant infections -infections on a host that is already infected. The differences between R E and R e are entirely due to the fact that R E ignores the fact that some people might already be infected. For malaria, R e could be computed in various ways, depending on the model, but at the steady state of any model with endemic transmission, R e = 1 by definition. We have not proposed any way of computing R e for this model.
• Local Adjusted Reproductive Numbers (or substituting R E for R C , as appropriate, to get Local Endemic Reproductive Numbers): D computes a measure of total HTC jointly spent among the patches, which describes how parasites are dispersed among patches by the human population strata (Eq 40).
G -The next generation matrix, |G| = (n + p) × (n + p), computed as a "type reproductive number." Note that most of the analysis is done with the diagonal blocks, R and Z, of the parasite's next generation matrix, G 2 .
J -Residency Membership Matrix, |J | = p × n. A matrix where the i, j th element is 1 if the j th population stratum resides in the i th patch. For example, J · H = P.
K -Mosquito dispersal matrix, |K| = p × p. The i, j th element of K is the fraction mosquitoes leaving the i th patch who end up in the j th patch.
By definition, diag(K) = 0 and j K i,j ≤ 1. Columns that do not sum to one have a net loss of mosquitoes associated with migration from the j th patch. R -The parasite's next generation matrix for human strata, |R| = n × n.
It is the upper on-diagonal block of G 2 , which describes transmission among strata in the parasite's next generation (Eq 45).
U -Egg Dispersal Matrix, |U (N , w ν )| = l × p. A matrix, derived from N and w ν describing how eggs laid by adults in a patch are allocated among the habitats (Eq 14).
V -Vectorial Capacity Matrix, |V| = p × p. A matrix that describes a patch-specific analogue of vectorial capacity. It describes the number of infective bites eventually arising in every other patch, on average, from all the mosquitoes blood feeding on a single person, on a single day, assuming all those mosquitoes became infected upon blood feeding (Eq 39).
X -Human Epidemiology State Space. A generic moniker for a state space defining for parasite infections and immunity in humans. The dynamics are given by a system of equations denoted dX /dt.
Y -Mosquito Infection State Space. A generic moniker for the state space defined for the component describing parasite infections of adult mosquitoes. The dynamics are given by a system of equations denoted dY/dt.
Z -The parasite's next generation matrix for patches, |Z| = p × p. It is the lower on-diagonal block of G 2 , which describes transmission among strata in the parasite's next generation among patches (Eq 45).

Upper Case Roman, A − Z
A -Ambient human population density, |A| = p. The density of humans present in a patch at a point in time (Eq 2) A r -Ambient resident human population density, |A r = (J ⊙ Θ) · H| = p describing the density of humans in each patch who reside in each patch, used as a diagnostic statistic for the time -spent matrix.
A n -Ambient non-resident, non-visitor human population density, |A − A r | = p describing the density of humans from the spatial domain who do not reside in the patch, used as a diagnostic statistic for the time -spent matrix.
B -Blood Host Availability (All Vertebrate Hosts), |B| = p. A measure of the of vertebrate hosts for blood feeding. It is used to compute the blood feeding rate, the human blood feeding fraction, and the emigration rate: D -Human Transmitting Capacity, |D| = n. The HTC is interpreted as an equivalent number of days (hence D) a person spends fully infectious [?]. It is computed by adding up days spent partially infections as a part of a day. To compute any spatial metrics involving reproductive success, D must be defined from the model X .
E -dEIR, described in the Abbreviations section above.
F -Functional Response F ν -Egg laying functional response, a function. F ν accepts a vector describing habitat availability, Q, and returns a vector of per-capita egg laying rates, ν. The function has a maximum rate ν x and a shape parameter s ν . See Eq 12.
F σ -Emigration functional response. A function that accepts vectors describing the availability of hosts, W , habitats, Q, and sugar sources, S and returns a vector of emigration rates, σ. It is assumed that mosquitoes move in search of resources, so the more available, the less the mosquito moves, but it will move in search of any one of three resources, so the rate is a sum of three terms, each with a maximum, σ * and a shape parameter, s * . See Eq 19. O -Other vertebrate host availability, |O| = p. O is a measure of the weighted, ambient population density of non -human vertebrate hosts term that is directly comparable to human availability, W .
P -The census population density, |P | = p. The number of people who live in a patch, called residents. Some patches might have zero residents.
Q -Aquatic habitat availability, |Q| = p. The availability of aquatic habitats in each patch is the patch-wise sum of the habitat search weights, w f . It is used to compute egg laying rates, ν, egg laying fractions, η, and mosquito emigration rates, σ. See Eq 11.
S -Sugar Availability, |S| = p. The availability of sugar affects mosquito movement. It is passed to the model as a vector of functions describing sugar availability over time.
W -Human availability, |W | = p. A measure of availability for blood feeding. It is a weighted measure of ambient population density, where the weights could take into account factors such as body size, ITN usage, or habitual use of personal protection. It is used in the computation of β, mosquito movement, including both the rate of emigration σ and the dispersal kernel, K.
W δ -Availability of Visitors, |W δ | = p. Availability of visitors from outside the spatial domain.
X -Infective Density of Infectious Humans, |X| = n. A dynamical quantity describing the effective infectious density of each stratum. This relationship must be defined in relation to the state space, X . Here, X = xH.
Y -Density of infected mosquitoes, |Y | = p. A variable describing the density of infected mosquitoes, defined by Eq 23.
Z -Infective Density of Infectious Mosquitoes, |Z| = p. A dynamical quantitity. Z is a variable describing the density of blood feeding, infective mosquitoes, defined by Eq 25.
Lower Case Roman, a − z b -Transmission efficiency, either b is a scalar or |b| = n. The fraction of infective bites that cause an infection.
c -Infectiousness, |c| = n. either b is a scalar or |b| = n. The fraction of infective bites that cause an infection.  f -Mosquito blood feeding rate, |f | = p. In each patch, the number of blood meals, per mosquito, per day (Eq 5).
g -Mosquito mortality rates, |f | = p. In each patch, the death rate of mosquitoes, per mosquito, per day.
h -Local Force of Infection, |h| = n. A term describing stratum specific hazard rate for local infection, the number of infections, per person, per day: The local FoI is computed from the dEIR h = f h (E). For the model family defined in this paper, h = bE.
k -Stratum infectiousness, |k| = n. A term describing the proportion of mosquitoes that become infected after blood feeding on a mosquito from each stratum (Eq 30).    w -Search weights w f -Blood feeding search weights, a vector of length n, are used to compute the availability of humans in each strata (See Eq 3). The search weights can be used to modify exposure to model usage of insecticide treated nets, differences in exposure related to body size or risky activities, use of personal protection devices, or other factors that modify blood feeding or exposure.
w ν -Egg laying search weights for aquatic habitats, a vector of length l. These search weights are used to compute the total availability of habitats, Q (see Eq 11). By using search weights to compute availability, it is possible to modify a baseline to consider various modes of control. Two possibilities are adding devices, such as oviposition traps, or fouling the water in a habitat.
x -Infectiousness, |x| = n. The probability a mosquito blood feeding on a human in each stratum becomes infected.
x 1 -True prevalence, |x 1 | = n. A term describing the true prevalence of P. falciparum infection in each population stratum (Eq 27) x 2 -Apparent prevalence, |x 2 | = n. A term describing the apparent prevalence of P. falciparum infection in each population stratum (Eq 29).
x δ -Infectiousness, visitors, |x| = n. The probability a mosquito blood feeding on a human visiting from outside the spatial domain becomes infected.
Θ -Time Spent Matrix, |Θ| = p × n. The j th column of Θ describes time spent by humans among patches.
The mosquito tracking matrix describes the location of cohorts of mosquitoes over time (Eq 37).
Ψ -Time at Risk Matrix, |Ψ| = p × n. The time at risk matrix modifies time spent by the mosquito circadian pattern, β -Biting distribution matrix, |β| = n × p. The matrix β describes the proportion of each bite occurring in each patch that is received by an individual in each stratum, defined by Eq 7. For example, (setting υ = 1), to compute the EIR for all the strata, we write β · f qZ. An important constraint is that the columns of diag(H) · β sum up to 1; β allocates bites to individuals in each stratum. To compute the FoI for mosquitoes (again, setting υ = 1), we write it in several different ways: γ -Net visitor infectivity, |γ| = p. In each patch, the fraction of infected mosquitoes that were infected by visitors.
δ -Travel FoI, |δ| = n. In each stratum, the net rate of infection from travel outside the spatial domain.
κ -Net infectiousness, |κ| = p. The proportion of mosquitoes in each patch who become infected after blood feeding on a human (Eq 10).
σ -Emigration rate, |σ| = p. The net immigration rate of mosquitoes from each patch is assumed to be a function of unsuccessful search for resources. See Eq 19.
ν -The population egg laying rate, |ν| = l. The egg laying rate by the adult mosquito population in each patch, per day (Eq 13). Eggs are distributed among the patches by availability (see η).
ξ -Circadian weight function, |ξ(d)| = 1. The function describes the relative blood feeding activity rates of mosquitoes. It is constrained such that 1 0 ξ(s)ds = 1 τ (t) -Extrinsic Incubation Period (EIP), |τ (t)| = 1. The number of days required for an infected mosquito infected at time t to become infective. We also define τ −1 (t), which is the time a mosquito becoming infective at time t became infected.
φ -Density independent mortality for larvae, |θ| = l. The death rate from all density -independent causes.
χ -A scalar. The number of eggs laid in a batch.
ψ -Maturation rate for immature mosquitoes, |ψ| = l. The rate immature mosquitoes mature. If ψ is constant then the time from egg to emergence is 1/ψ.